% Darzi et al. (2022), Short-Term Bayesian ETAS Spatiotemporal Forecasting of the Ölfus 2008 Earthquake Sequence in Iceland
% submitted to Journal of Geophysical Research: Solid Earth
% Atefe Darzi (atefe@hi.is)

% To regenrate Figure 7a: 
% Daily aftershock seismicity forecasting map for earthquakes with M_w≥2 over the aftershock region 
% following the 29 May 2008 mainshock at 15:45 UTC (Mw6.3). 
% The forecast is issued at 00:00 UTC on 30 May 2005. 

%% INPUTS and DATA DESCRIPTION
% spatialN_mgrM.mat file includes the number of forecasted events with M>=m
% over each spatial grid cell across the aftershock zone corresponding to each MCMC
% samples of ETAS model parameters (sampled posteriors)
% num_grid: Number of grid cells in aftershock zone: 7301
% 198: Number of MCMC samples obtained from Bayesian parameter estimation technqiue

%% OUTPUTS
% plot_ErrorBar: The median number of forecasted events shown in the gray box 
% with the posterior distribution: 16th and 84th percentiles (blue numbers) and 2nd and 98th percentiles (black numbers). 
% plot_ForeSeis: aftershock seismicity forecasting map (median estimates) for earthquakes with M2+ 

%% Seismicity Forecasting Map
% 1998-2008 EQ aftershock zone 
lonMin = -21.6;
lonMax = -20.8;
latMin = 63.74;
latMax = 64.23;
% spatial grid cells
ratio = 0.0151/0.00664; % for Iceland
deltaGrid_Y =  0.01; 
deltaGrid_X =  ratio*deltaGrid_Y; 
Xgrid = lonMin:deltaGrid_X:lonMax;
Ygrid = latMin:deltaGrid_Y:latMax;
   Xcgrid = Xgrid(1:end-1)+deltaGrid_X/2; % long
   Ycgrid = Ygrid(1:end-1)+deltaGrid_Y/2; % lat
   num_grid = length(Xcgrid)*length(Ycgrid);

% lower cut-off magnitude
Ml_cat = 2.0; 
% magnitude thresholds 
vec_m = [Ml_cat 3.5:1:6.5];

load('spatialN_mgrM.mat') % sampleN_mgrM

% Calculate 2nd, 16th, 50th, 84th and 98th percentile of forecasted number of events for M> m 
N_robust_p50_xy  = zeros(length(Ycgrid),length(Xcgrid),length(vec_m));

for k=1 :length(vec_m)
% for the whole aftershock zone
    [N_robust_mean(k),N_robust_p50(k),N_robust_p16(k),N_robust_p84(k),N_robust_p02(k),N_robust_p98(k)] = ordered_statistic(sum(spatialN_mgrM(:,:,k))); 

% over each grid cell for spatial distribution of forecasts
    [N_robust_mean_grid,N_robust_p50_grid,N_robust_p16_grid,N_robust_p84_grid,N_robust_p02_grid,N_robust_p98_grid] = ordered_statistic(spatialN_mgrM(:,:,k));
    N_robust_p50_xy(:,:,k)  = (reshape(N_robust_p50_grid ,length(Xcgrid),length(Ycgrid)))';
end

P_mean(:) = 1-exp(-N_robust_mean(:));

%  RUN the codes below 
plot_ForeSeis
plot_ErrorBar






    
    